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ABSTRACT 


The challenge m prediction of compressible flow phenomena is to resolve as accurately as 
possible the sudden discontinuities of the flow vanables associated with shocks Over the last 
three decades, extensive research has been earned out in this direction Different algonthms 
have been developed to solve Euler and Navier-Stokes equations for compressible flows The 
algonthm developed so far can broadly be classified into two categones One is of central 
difference type with explicitly imposed artificial dissipation such as MacCormak and Jemson 
type algonthms The other is the upwind type namely Steger and Warming, Van Leer 
splitting where the discretization is made followmg the physics of flow and the direction m 
which the information is propagated It has been thought that the upwind scheme would 
perform relatively better as it is based on the physics of the flow The flux vector type of 
splitting IS one that is simple to implement but leads to an maccurate solution because of its 
high dissipation properties This high dissipation can be eliminated by higher order 
discretization The Flux difference splitting has shown improved accuracy such as Roe's Flux 
Difference Splitting scheme However, the flux difference schemes are complex to 
unplement and takes more computational time Recently a scheme named as "Advective 
Upstream Splitting Method (AUSM)" has been developed combining the advantages of both 
of the above two schemes The AUSM algonthm has been found relatively simple for 
codification and also minimum as accurate as Roe's Flux Difference Splitting Therefore m 
the present work an attempt is made to develop a flow solver usmg the AUSM algonthm m 
the fimte volume formulation The Navier Stokes equations are mtegrated m the physical 
plane so that the truncation error associated with the transformation of vanables from the 
physical to computational plane can be eliminated To split the fluxes at the cell faces 
following the algonthm, the vanables from the cell center are extrapolated by piecewise 
fimte element reconstruction method The logic is extended to formulate the higher order 
discretization To judge the performance of the flow solver, different test cases such as flow 
over a bump with subsomc and supersomc conditions, both mviscid and viscous flows, flow 
over a flat plate with shock boundary layer mteractions are computed The predicted results 
illustrate the flow solver developed is accurate enough m resolvmg the boundary layer 
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growth, formation of shock associated with compressible flows with an accuracy acceptable 
for engmeenng applications 

Another objective of the present thesis work is to find the applicability of the flow solver to 
predict the complex nature of the turbulent flow by Large Eddy Simulation (LES) The huge 
computational resource required in performing a direct numencal simulation with moderate 
Reynolds number applications and over complex engmeenng geometnes makes it a distant 
approach The practical engmeenng problem can now be tackled with LES where the Larger 
eddies are resolved directly and the smaller eddies, which are universal m nature, are 
modeled by a suitable Subgnd Scale Model Thus we have incorporated LES m the present 
solver to resolve the pattern of flow over a backward facing step This is to judge the 
suitability of the flow solver in LES calculation However for the limited computational 
facility available the analysis is restncted only in two-dimension It has been seen that the 
solver IS able to predict the vortex structure generated downstream of the step Further the 
averaged velocity field has been found to match reasonably well with the expenmental data 
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CHAPTER -1 


INTRODUCTION 


The flow fields m engmeenng applications are associated with turbulence, mixing, separation 
and transition An appropnate flow solver together with an accurate turbulence model is 
always necessary to predict these kinds of flow and heat transfer phenomena Even with the 
sophisticated turbulence model and solution algonthms, it is very difficult to predict properly 
the sharp changes of flow vanables across the shock in compressible flow, the thin shear 
layer and the flow separation The conventional turb> fence models are good for certain 
engineering applications, but they are not universal in ature One way is to use the Direct 
Numencal Simulation, which is restncted to only lou Reynolds number applications and on 
simple geometnes because of its high computational cost Here lies the importance of the 
Large Eddy Simulation (LES) where flow properties can be evaluated by solving the filtered 
Navier-Stokes equations in finer grids and modeling the turbulent stresses by a suitable 
Subgnd Scale Model However, an appropnate flow sciver is needed for LES calculation. 
The flow solver should be such that it can resolve the complex flow features due to shock- 
boundary layer interactions and sudden discontinuities that occur in compressible flows with 
minimal dissipation level A well-designed flow solver preferably based on an upwind 
scheme can atleast minimizes the problem The widely used certain central difference 
algonthms with explicitly imposed artificial dissipatici (like Lax-Wendroff, MacCormak, 
Jameson algonthms) may not be suitable for LES At the same time the dissipation level m 
some Flux Vector Splitting algonthms may corrupt the solution near the some pomt or withm 
the VISCOUS sublayer On the other hand, though the Roe's flux difference splittmg (1982) is 
reasonably accurate, it is quite difficult for codification and to extend for a three-dimensional 
problem Recently a scheme named as "Advective Upstream Splittmg Method (AUSM)" 
(1993) has been developed which has shown to overcome most of these difficulties This 
scheme is relatively easy to implement and needs a lower amount of computation. The 
AUSM seems to be promismg even over the Roe's flux-difference-splittmg scheme 
Therefore an attempt is made to develop a two-dimensional flow solver m fmite volume 
formulation by using the AUSM algonthm 
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dimensional compressible Navier-Stokes flow solver m the finite volume formulation 
based on the AUSM algonthm The attractive features of it are that equations are 
integrated in the physical plate so that any complex geometry can be handled properly, 
and less truncation error would be resulted ITie present work is based on the flow solver 
developed by Sarkar (2000), where the Navier-Stokes equations are solved in the physical 
plane based on the four-stage Runge-Kutta algonthm in finite volume formulation. The 
Jameson algonthm used to split the convective fluxes of the Navier-Stokes Equations has 
been changed by AUSM algonthm A piecewise finite element reconstruction procedure 
has been used to extrapolate the flow vanables at the cell faces so that the solver would 
be at least second-order accurate in space A local transformation procedure has been 
used to evaluate the viscous stresses in the Navier-Stokes equations The local time 
stepping and residual smoothing have been used to accelerate steady state solution of 
some of the test cases 

Vanous test cases have been performed to assess the performance of the flow solver for a 
wide range of flow conditions It has been seen that the flow solver is able to resolve the 
shock wave and the separated flow with reasonable accuracy The objective of a new flow 
solver was to study the applicability the solver for LES calculation Thus the solver is 
used along with a Subgnd Scale Model to predict the flow over a two-dimensional 
backward facing step with 8 9 expansion ratio The analysis is done for a very low firee- 
stream Mach number (M„=0 128), for which expenmental data was available As it is 
found that in the near wall region the Smagonnsky model is highly dissipative in nature, a 
damping function is used in the subgnd scale stress term to minimize the near wall 
dissipation The predicted results have been found to mateh reasonably well with the 
expenmental data. 

The thesis consists of four chapters m addition to the first chapter, Introduction, as 
mentioned below 

A cntical review of the existing literature related to different upwind schemes and LES 
models is presented and documented in Chapter 2 Chapter 3 presents the governing 
equations and the solution technique in finite volume formulation The thorough 
explanation of the Smagonnsky model and its modification to capture near wall 
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analysis of the applicability of the upwind solver is presented The second part of the 
chapter deals with the computation of flow over a backward facing step by LES 
Conclusions based on the present study and some thoughts on the direction of future 
research work in this field are provided in Chapter 5. 
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CHAPTER - 2 
Literature Review 


The objective of the work is to develop a flow solver based on the charactenstic 
information of flow so as to use the solver for LES and DNS type calculations Thus, m 
the present chapter a review of the literature for various upwind schemes is presented 
Some studies on LES computation and subgrid scale models are also discussed 

2.1 Upwind schemes: 

The first upwind scheme was developed by Courant etal (1952) based on the 
charactenstic form of the hyperbolic equations, where the discretization is based on the 
sign of the eigenvalue The central difference of the hyperbolic equation is unstable near 
the sonic transition point It can be shown that one-sided difference leads to the stability 
of the hyperbolic equation, when the difference is done based on the charactenstic 
eigenvalue of the equation The scheme of Godunov (1959) comes from the local exact 
solution of the Euler equations, which is stnctly conservative m nature. Here the solution 
was considered as piecewise constant over each mesh cell at a fixed time and the 
evaluation of the flow at the next time step results from the wave interactions at the 
boundanes between the adjacent cells The solutions to the Riemann problem at each cell 
interface produces a perturbaton of the fluid state resulting from the propagating waves 
over the time interval At 

The scheme developed by Warming and Beam (1976) is a second-order upwind scheme 
The upwinding of the convective fluxes is performed by modification of the corrector 
term of the MacCormack scheme One-sided difference is used instead of the central 
differencing whenever the local charactenstic speed was found to be of the same sign It 
was found that switching from fully MacCormack to the upwind scheme loses the stnet 
conservation A spatially transmission operator was denved for switching from the former 
to the latter to make it fully conservative in nature This is essential for captunng of the 
shock waves in the flow domain However it was found that the scheme has supenor 
dissipative and dispersive properties than those of the central schemes. 
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Monetti (1979) developed an explicit first-order upwind algonthm by recasfing the Euler 
equation to the charactenstic (Riemann) type vanables This algonthm was known as 
explicit X scheme A two-step MacCormack algonthm was used for the time stepping 
The vanables were discretised by two point windward difference for negative eigenvalue 
and three point windward difference for the positive eigenvalue at the predictor level, 
along with three point windward difference for the negative eigenvalues and two point 
windward difference for the positive eigenvalues at the corrector level Though it was 
found that the scheme was well suited for the transitional flows, it was not found to be 
conservative in nature The scheme was found to be complicated for the two- and three- 
dimensional flows and geometry dependent which was not relevant from the physical 
point of view 

The work of Osher and Enquivist (1980) descnbes the conservative upwind Riemann 
solvers, which are monotomc in nature They found it suitable for transonic flows with 
small disturbances where the compression and the rarefaction waves are used to 
approximate the shock waves The work was further modified by Chakravorty and Osher 
(1983) to analyze a method for treating computational and physical boundanes based on 
solving initial boundary value Riemann problem The one dimensional flow, quasi-one 
dimensional Laval nozzle flow, conical supersonic flow past an aerofoil and Mach 8 
supersonic flow past a circular cylinder were simulated by this algorithm and were 
presented by them in the same year (1983) The results were found to be quite agreeable 
with the available data 

A remarkable property of the flux vectors of the inviscid gas-dynamic equations is that 
they are homogeneous functions of degree one Using that properly Steger and Warming 
(1980) denved the flux vector splitting algonthm of mviscid gas-dynamic equations by 
one-sided differencing of the fluxes depending on their charactenstic eigenvalues Both 
explicit and implicit finite difference algonthras were denved In the explicit algonthm, 
MacCormack algonthm was modified in such a way that tiie forward difference is applied 
to the predictor step and backward difference in the corrector step and only one-sided 
differencing was allowed to make it a completely upwmd algonthm Thus the stability 
bound was increased to 2 for one-dimensional case Several numencal experiments on 
one and two-dimensions were performed and found to be agreeable with the actual 
results However it was found that the split fluxes are not always continuous and shows 



non-differentiabihty at some transition points This happens because of the eigenvalue of 
the charactenstic equation changes its sign leading to the numencal oscillations near the 
transition or at the sonic points These oscillations can be reduced by introducing a small 
dissipation parameter in calculating the eigenvalues of the charactenstic equations 
(Burning P G and Steger J L (1982)) The more efficient approach is to split the fluxes m 
such a way that the slopes of the flux components tends to zero at the transition point, so 
that no discontinuity would be maintained at the some point as suggested by Van Leer 
(1982) In such cases the fluxes and the associated Jacobians are made continuous 
functions of the Mach number and expressed as the polynomial of the lowest possible 
order The scheme was found to be more robust in captunng the strong shock 
discontinuity as compaied to Steger and Warming flux vector splitting scheme 

Dadone and Napolitano (1983) modified the explicit lambda model to an implicit one for 
one-dimensional and two-dimensional flows The compressible conservative equations 
were expressed in the form of charactenstic (Riemann) vanables The old time level 
denvatives were discretised by three-point second order accurate windward difference 
and unknown denvatives were discretised by two-point first order accurate windward 
difference Therefore the resultant discretised equations were converted to a block 
tndiagonal system of equations for one-dimensional flow However to retain the block 
tndiagonal form for two-dimensional flows they solved the system of equations by ADI 
scheme The model was further modified to a three-dimensional system of equations for a 
general orthogonal curvilinear coordinate system However the method appeared 
computationally expensive 

The work of M Ben Arzi et al (1984) describes the solution of generalized Riemann 
initial value problem with the help of second order accurate Godunov's method They 
analytically denved the time denvatives bofli in Lagrangian and Eulenan frame of 
references even for the flow with the jump discontinuities The scheme is found to be 
conservative and second order accurate in both tune and space. 

One of the Godunov types of scheme named as “Flux Difference Splitting” developed by 
Roe (1981) was an approximate Riemann solver. The fluxes in the Riemaim equation 
were decomposed following the charactenstic speed of the wave propagation The 
upwmd dissipation term was calculated by usmg the local Jacobian matnx of the flux 
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terms The vanables used to compute this matnx were denved by an average weighted by 
the square root of the local densities across two sides of a cell face The scheme was 
found to be first order accurate in space and stnctly conservative in nature Another good 
feature of the scheme was that it was able to capture the discontinuity and shock more 
accurately than flux vector splitting algonthms 

The interesting feature of the Steger and Warming flux-vector splitting is its simplicity m 
implementation and that of Roe flux-difference splitting is its accuracy These properties 
were combined by Liou and Steffen (1993) to get a new type of scheme named as 
“Advective Upstream Splitting Method (AUSM)” The convective and the pressure terms 
were separately discretised in this method The upwmdmg of the convective terms was 
made depending on the cell interface velocity, which was obtained by splitting the Mach 
number at the cell interface following the Van-Leer splitting The pressure was weighted 
by using the polynomial expression of the charactenstic speeds (M±l) The scheme is 
relatively easier to implement and also needs lower amount of computation time than 
FDS scheme Further the tested results indicated that the scheme could be implemented 
over any complex geometry with accuracy as high as Roe’s FDS scheme The dissipation 
level and the entropy generation within the viscous layer in minimal even compared with 
the Roe’s scheme and thus it resolved correctly the nose region a supersonic blunt body 
It has been felt the scheme would be a suitable for LES/DNS type calculation for complex 
geometnes 

2.2 LES computations: 

The DNS results appeared so far m the open literature are restricted to simple geometnes 
and for low Reynolds number flows due to high computational time An alternative is to 
perform LES calculation where the large scale eddies are captured and the smaller scale 
eddies which are more universal m nature are modeled A bnef discussion on LES and 
subgnd scale (SGS) models are presented The governing equations of turbulent flows are 
given in LES by filtermg the Navier-Stokes equations The turbulence can be thought of 
superposition of wakes with different wave numbers The large-scale motions generated 
by the larger wakes are flow dependent where the small-scale wakes are found to be 
universal m nature To eliminate the subgnd scale fluctuations it is necessary to apply a 
filtenng operation on Navier-Stokes equations that essentially cuts off the wavelength 
smaller than Kc (the charactenstic wave number defined to be ir/A, from the Nyquist 
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critenon, where A is the filter width) In that sense, for a large eddy simulation at a 
practical Reynolds number should have a filter cut off at a wave number m the inertial 
range SGS models should model the range at which energy decays exponentially 

2 2 1 Filtenng approach 

Different filtenng operations are required for different levels of turbulence Germane 
(1986), Mom and Kim (1982), Germane (1992) tested many filter previously The 
Gaussian filter was found to be most suitable because it has the same expressions as in 
physical as well as m Founer space Erlebacher etal (1992) and others made the similar 
conclusions But for the wall-bounded flows vanable length filters should be adopted 
because of the vanation of turbulent length scale perpendicular to the wall as suggested 
by Mom etal (1982) Several other filters were tested among which top hat filter is 
widely used because of simplest expression m Founer space 

2 2 2 Modellmg of SGS terms 

A basic ingredient that charactenzes the subgnd models used m LES is that their 
charactenstic length-scale is given by the gnd size As a consequence simple algebraic 
models or one-equation turbulence models were adopted The simplest one is 
Smagonnsky model (1963) The model was adopted by correlating the anisotropic 
deviatonc part of the stress to the kinematic diffiisivity It was found that the 
Smagonnsky model of eddy viscosity has a one to one correlation with the local strain 
rate tensor and it is also found to be independent of the filtering operation The particular 
type of averaging was reflected over C* (the Smagonnsky constant) and A For isotropic 
turbulence the value of Cs was analytically determined by Lilly m 1987 From his 
calculation Cs was found to be 0 18 But m general to avoid excessive subgnd scale 
dissipation Cs is limited to 0 1 This makes the model to behave reasonably well for firee- 
shear flows It was also seen firom the work of Deodorff (1970) and Schumann (1975) that 
this model is able to predict the charactenstics of turbulent channel flow However, it was 
found that m the near-wall region where the viscous stress is predominant the 
Smagonnsky model is too much dissipative m nature So the Smagormsky model failed to 
detect the viscous sublayer and the buffer layer for a wall bounded flow In such cases the 
wall model should be adopted and some dampmg functons should be used to make the 
eddy-viscosity negligibly small very near the wall The following modification of the 
Smagonnsky constant can be as, Cs=0 l[L0-exp {-(y‘^/25)}^]*^^ where y^ is the non- 
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dimensional distance away from the wall Bui (1999) suggested the same approach to 
simulate compressible flow inside a square duct In 1992 Amal and Fnednch used a 
different version of this model by replacing CjA with mm (CsA, Im), where U is the 
mixing length in the near-wall region denved from Nikuradse (1933) The Sraagonnsky 
constant coupled with the damping function was also used by M Breuer and W Rodi 
(1994) to analyze flow over a 180° bent duct, and B J Boersma et al (1994) for the flow 
simulation inside an electronic flow-meter Similarly others adopted the same procedure 
like T G Thomas etal (1994) where CjA was replaced by min(CsA,k Imr(lw^)) and 
r (lw^)=l 0-exp(-lwVA'^), and M Manhart etal (1994) where CjA is replaced by min(k 
Xn and Cs A ) and Xn was the normal distance for simulation of flow over a cylinder It is 
also possible to use different values of Cs at different time as suggested by A Andren 
et al (1994) Thus, Cs is to some extent flow dependent 

2 2 3 Scale-similantv and mixed Models 

The eddy viscosity closure assumed a one to one correlation between the subgnd scale 
and the large-scale strain tensor However m actual cases the correlation between these 
two quantities was found to be very less This lack of correlation made Bardina et al 
(1980) to propose an alternative subgnd scale model called scale-similanty model This 
model is based on the double filtenng of the flow vanables The basic idea originates 
from the interaction between the smallest scale eddy of the resolved scale and the largest 
scale eddy in the subgnd scale that causes the energy back scatter from the subgnd scale 
to the resolved scale The inclusion of the Leonard and Cross-term in the subgnd scale 
tensor made the model capable m the prediction of subgnd scale backscatter of the kinetic 
energy Indeed the good correlation with the DNS data made the model highly important, 
which was further modified to the dynamic model 


2 2 4 Dynamic SGS modelling 

The basic idea of the dynamic SGS model, developed by Germane etal (1991) was to 
improve the Smagonnsky model by setting Cs different values for different mesh nodes at 
each time step, which attempts to self-adjust the dissipation by subgnd scale stresses 
This model was developed basically by double filtenng of the Navier-Stokes equations 
At first a gnd level filter was used to resolve the fluctuatmg velocity field Next a test 
level filter that was generally larger in magnitude was applied on the filtered velocity 
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field The excess of fluctuating resolved velocity field and filtered resolved velocity field 
obtained firom test filter gives the information of the value of the constant 

Cdyn (x, y, z)= , where Lm is the Leonard term denoted as 

2MkiSki 

Cl Cl 

Lk!= W,w, -U,M, 

Mki = ^ 1 5 1 iS;;; - Z? 1 5 I S'y , where Ski is the local strain rate and | S j = (2 

However the model was proved to be mathematically inconsistent by Ghosal etal (1995), 
because C was taken to be independent of the test filter width From the test using channel 
flow data from direct numencal simulations, Germano etal (1991) shown that the 
denominator of the expression could locally vanish or become sufficiently small, thereby 
leading to computational mstability. Later Lilly (1992) adopted a least square procedure to 
mmimize the error so that the expression of C can be given as 

2MkiMki 

Sometimes this dynamic constant causes the instability of the simulation as it vanes in a wide 
range over the whole flow domain producing a strong gradient of Ut. It is also possible that 
at a particular gnd pomt the constant remains negative for a larger penod of tune It leads to 
subgnd scale energy transfer to the resolved scale causing the excessive backscatter of 
energy To avoid this some temporal and spatial averaging may be done along the 
homogeneous direction The simplest way to modify the constant suggested by P. Sagut et al 
(1994) IS to truncate the dynamic constant by the followmg constraints 

1 — +ut>0 

Re 

2 Cdyn — Ctnax 

The maximum value of the dynamic constant is taken as the square of Smagormsky constant 
This model can be able to predict the near wall boundary layer situations In another way to 
restnct the value of C to a time mdependent function and to allow variations with low 
frequencies, a low pass filtering method was chosen (M Breuer and W Rodi (1994)) C can 
be given as 

Cfi,tered"^‘=(l O-S) C" + 8 


1 A 



With value of e ~10‘^ All frequencies of oscillation damps out with low frequencies of 
oscillation remams Akselvoll and Mom (1993) also adopted the method for backward 
facing step, which has still one homogeneous direction 

The mathematical inconsistency presented in Germano’s model was removed by the 
newly developed constrained dynamic localization model as proposed by Ghosal et al 
(1995), which can be imposed on flow without any homogeneous direction The constant 
C > 0 was m this model m order to denve an integral equation for C However, this model 
shows no backscatter, which is essential for simulation of transitional flows Another 
model was proposed by the same author (known as k-equation locali 2 ation model) can be 
able to transfer subgnd scale energy to the resolved scale by the inclusion of the transport 
equation for subgnd scale kinetic energy However these two models are still not widely 
used for inclusion of two more integral equations and one transport equation 

2 2 5 SGS for compressible flows 

Yoshizawa (1986) suggested the first compressible SGS model for weekly compressible 
shear flow with the aid of multi-scale direct interaction approximation The model had 
included Leonard and cross terms that are generally neglected m the Smagonnsky model 
In fact the effect of back scatter due to cross-term is not clear for the Smagonnsky model 
This IS only applicable for the cross term because of the use of an asymptotic expansion 
about an incompressible state The model was further modified to the compressible flow 
calculation by Erhebacher etal (1992), Vreman etal (1994), Squnes etal (1991) and 
Mom etal (1991) The model presented by Erhebacher etal contains the crossed and 
Leonard terms both for modelling of the subgnd scale part of the momentum and energy 
equations, which was otherwise neglected by Yoshizawa’s model Though the isotropic 
part of the turbulent stress tensor was modeled it was found that the model correlates 
poorly with the DNS of isotropic turbulence presented by Speziale et al (1988) Bui 
(1999) confirmed the same result with inclusion of this term m the LES model However, 
with the low turbulent Mach number, which is generally encountered m the practical 
cases, the isotropic part of the turbulence is an order of magnitude lower than the 
thermodynamic pressure So, it is generally neglected m the turbulent flow simulation of 
low Mach number flows This was also supported by Squines et al. (1991) and Vreman 
et al (1994) This model was further improved Mom et al (1991), Squnes et al (1991), 
where the Smagonnsky constant was determined dynamically m space and time so that 
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the model can be made flow independent However the models are restricted only within 
the low Mach number range and simulation of high Reynolds number and high Mach 
number flows are still not found in open literature 

2.3 Conclusions: 

The present review of literature reveals that a significant progress has been made in the 
development of upwind schemes and different t3Tjes of LES models Applications of 
these models on relatively simple geometnes illustrate that the models are successful to 
predict the largest structure of turbulence with high accuracy However the simulation of 
high Reynolds number flow over the complicated geometnes is still a formidable task and 
there is further scope of improvement m this area 
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CHAPTER -3 


GOVERNING EQUATIONS AND NUMERICAL PROCEDURE 

3.1 Introduction 

The mam objective of the present work is to develop a two-dimensional Navier-Stokes 
solver by using the AUSM algonthm Also the present study explores the ability of the 
Smagonnsky LES model to simulate the flow over a backward facing step Thus, m this 
chapter, both the Euler and the Navier Stokes equations, the LES model and the 
numencal technique in discretisation of the governing equations are presented 

3.2 Governing equations and LES model 

In this chapter, the Euler equations and the Navier-Stokes equations for the laminar and 
turbulent flows are presented All these equations are presented m the conservative forms 
and solved m physical plane based on finite volume algonthm For the turbulent flow, the 
governing equations can be obtained by filtenng (local volume averaging) the Navier- 
Stokes equations Further to reduce the equation to a simpler form a density weighted 
averaging (Favre filtenng) has been done The density weighted averaging leads to a 
simpler form of the governing equations with terms that are more amenable to physical 
interpretation The Subgnd Scale Stresses are modeled by the Smagonnsky model It 
should be noted that density-weighted space averaging is used for velocity components 
and temperature, whereas the space averaging is used for pressure and density 


3 2 1 Flow equations 

The time-dependent, two-dimensional, compressible Navier-Stokes equations for laminar 
flows can be expressed m the conservative form as, 


dt dx dy 


= H 


(3 1) 


where. 


fp ] 


"pu ^ 

pu 


pu^+p-txx 

pv 

, F = 

pUV-Txy 

U; 


^(E + p-txx)U-txyV + qx, 
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r 


pv 

PUV-Tyx 
2 



^0" 


0 

and H= 

0 


loj 


(3 2) 


pv +p-Tyy 

(E + p-Tyy)V-TyxU+qj ^ 

are state vector, convective flux vectors m x and y directions and source terms vector 
respectively The stress tensor and the heat flux vector are given in Cartesian coordinates 
as. 


Ty — |i.| 


au, 5u, 


^axj 


5x, 


_2g auk 


3 yaxk 


and q, 


i^Cp ex 

Pr axj 


(3 3) 


In the present study, the laminar viscosity coefficient p. is assumed to be a function of 
temperature only, and is evaluated following the Sutherland’s law which is given by, 
v3/2 


p = C 


(T)- 


-, , where Ci and C 2 are constants for a given gas For air at moderate 


T + C 2 

temperatures, Ci= 1 458x10'® kg/(ms V°K) and Ci= 1 10 4°K The perfect gas equation of 
state p = pRT , is considered to be applicable 


For an mviscid flow, Eq 3 1 remains the same where the elements of the column vectors 
are simplified Thus, omitting the viscous terms from the Navier-Stokes equations, the 
Euler equations in the conservative form are, 


dU dF dG 
dt dx dy 


(3 4) 


where. 




^pu '' 


fpv ) 

pu 

F = 

pu^+p 

and G = 

puv 

2 

pv 


puv 


pv 4-p 

U; 


v(e+pK 


<(E + p)v, 


(3 5) 


3 2 2 Governing equations for LES 

The governing equations for compressible turbulent flows can be obtained by filtenng the 
Navier-Stokes equations The equations are similar to that of the laminar flow equations, 
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but there is additional subgnd scale stress terms m momentum and energy equations 
They can be wntten as 

^ + ^ + ^ = 0 
3t 5x 5y 


rp ] 


1 






pu +P-T3,^ 

, F = 


pv 


PUV - X,^y 



^(E + p-Xxx)u-x^yV + q^^ 


U = 


pv 

puv- V 

+P-Tyy 

(E + p-Tyy)V-XyxU+qy 

The stress tensor and the heat flux vector can be reduced to, 

l^Cp af 


and 


(3 6) 


-r,, = P- 


^5u, 5u,^ 

L + L 

9x, 


3 ^ 5x]j^ 


+ Tt,j and q, 


Pr 5X] 


+ qti 


(3 7) 


The bar in the equations (3 6) and (3 7) denotes a flow quantity, defined as 

f = jG(x-x')f(x')dx', where G is a special filter and the integral is over the flow 
D 

domain D The tilde quantities in the flow equations can be obtained by Favre-filtenng 

7 pf 
given as 1 = -^ 

P 

The turbulent stress quantities Xt,j and % have been modeled by a Subgnd Scale Model 
(SGS) model as used by Bui (1999) It has been discussed m detail in the next section 


3 2 3 Subgnd Scale Model 

The subgnd stress term in the momentum equation is 
where, 

x,j =(u,Uj-u7uj)f 

This term can be modeled by 

= 2CpA^ I S I (S,j - j W„) (3 9) 

where. 
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q^=^t. 


IS the trace of the SGS Reynolds stress tensor The filtered velocity gradient tensor is 


5u, 

V^J 


+ - 


duj 

9x, 


and 

|SK2SyS„)'^' 

In Eq 3 9, C IS a constant to be determined according to the particular SGS model used 
For Smagonnsky SGS model, a value of C=0 1 is chosen which is the square of the 
Smagonnsky constant Cs 


For the wall bounded turbulent flow, the Smagorinsky model was found to be too much 
dissipative near the wall So this Smagonnsky model has to be modified near the wall to 
capture the viscous sublayer there A typical way is to vary the Smagonnsky constant 
dynamically throughout the flow domain But this process was found to be too much 
computationally inefficient even for simpler flows A simpler way is to multiply the 
Smagonnsky constant with the Van-Direst damping function to make the eddy viscosity 
negligibly small near the wall The resultant constant is 


C=0 1(1 0-exp 


f 


3^ 


y 



25 


1 

k J 

) 


) 


(3 10) 


Where d is the normal distance from the wall in wall units, defined as 

and the fiiction velocity is defined as 


Ut= 


^ W 


After the step the normal distance is calculated firom two wall directions and is reduced to 
be 


d= (3 1.) 

The filter width A is chosen to be the gnd spacing size, which is defined as Va , where A 
IS the area of the cell 
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From a pnon test done by Vreman et al it was found that the inclusion of the SGS 
isotropic stress tensor leads to instability of the system of equations So in our study we 
have neglected the isotropic part of the turbulence With the term omitted, the SGS 
stresses included in the momentum equations by simply replacing the laminar viscosity 

coefficient pgff = p. + pscg where, psQs = CpA^ | S | , and in the energy equation it is 


included as qefr=q+qsGs, where qsGs = - 


l^SGsCp 5f 
Pq 9x, 


3.3 Numerical scheme 

The unsteady, two-dimensional Euler and Navier-Stokes equations are solved in the in 
finite volume by using the four-stage Runge-Kutta time-stepping scheme, in conjunction 
with different accelerating techniques 


The computational domain is divided into quadnlateral finite volume cells, Fig 3 1 For 
each cell, the unsteady two-dimensional Navier-Stokes equations Eqs 3.1 can be 
represented m the integral form as. 


at 


JUdQ+ jQdS=0 

a s 


(3 12) 


where, Q. denotes the volume of the control volume being considered, dS is the surface 
normal, Q = (F,G) is the flux vector and U is the state vector 
A simple discretised form of the Eq 3 12, for a discrete cell (i, j) Fig. 3 1, is given by, 

|-(a,,jU,,j) + <f(Fdy-Gdx) = 0 (3.12a) 

^ S 


£ 

dt 


(^,jUx,j) + Q(U),,j 


= 0 


(3 12b) 


where, Q,j is the volume of the discrete cell and Q (U)y are finite volume approximations 

for the net convective and diffusive fluxes out the discrete cell The mdexmg scheme and 
the integration direction of the line integral is shown in the Fig 3 1 As an example for 
the east face of the cell the values of Axi and Ayj are given as 

Axi = Xj+ij+i - x,+ij 
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^yi Yi+ij+i “ Yi+ij 

and Fk and Gk are x and y direction flux averages on the kth cell face given as, 

Fk=[pu, puu+p, pvu, (E+p)u 
Gk =[pv, pvu , pvv+p, (E+p)v 

The integral QdS m Eq 3 12 is thus taken care of by the combination of the convective 
and diffusive fluxes out of the discrete cell For a cell-centered scheme, U,j is assumed to 
be at the center of the cell The convective term is computed my modifying the 
“Advective Upstream Splitting Method”, developed by Liou and Steffen (1993) for the 
arbitrary finite volume cells and making it upto the second order accurate by tihe method 
developed by Coiner (1994) For the diffusive flux, the calculation procedure of the 
viscous denvatives is discussed m a later section 

3 3 1 Calculations of the convective fluxes 

For the calculation of the advective fluxes the solution procedures are reconstruction and 
flux computation The pnncipal aim of the reconstruction procedure is to make higher 
order upwindmg of the convective terms For this reconstrachon procedure the pnmitive 
vanables (p, u, v, p) are calculated at the cell center and extrapolated at the cell faces to 
discretize the convective terms For the second order upwindmg each of the four pnmitive 
vanables p, u, v, p are assumed to vary linearly within the finite volume as 

(j) x= $ +Ux (x- X )+Uy (y- y ) (3 1 3a) 

where, <}) will be any of the above four vanables For the third order upwindmg process 
these vanables are calculated as 

(j) x=^ +Ux(x- X )+Uy(y- y )+Uxx(x- x )^/2+Uyy (y- y )^/2+Uxy(x- x )(y- y ) (3 . 1 3b) 

Following these extrapolation formulae the values of the vanables can be extrapolated to 
the cell faces The bars in equation (3 13a and 3 13b) denote the cell averaged value and 
Ux, Uy, Uxx, Uyy uud Uxy dcnote the gradients withm the cells 

The gradients Ux, Uy are computed by the method used by Comer (1994) The value of 
Ux, Uy in the target cell is computed by the least square procedure that m immizes the sum 
of the square of the differences between the values of reconstmction polynomial of tihe 
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neighboimng cells and that of the target cell For a two dimensional confrol volume the 
minimization procedure can be wntten as 


’aibi‘ 

Ux' 


'ci ' 

.^1^2. 

Uyj 


_C2_ 


aj = -xo)^ 

1=1 

4 

^ 2 = -Yof 

1=1 

4 

= ~yo) 

1=1 

E(Ui“Uo)(Xi -Xo) 

1=1 

c2=i;(lJ,-Uo)(y,-yo) 

1=1 

where, i=l-4 denotes the neighbounng cells and i=0 denotes the target cell These 
gradients are used to extrapolate the values of the dependent vanables on the cell faces 
The values of Uxx, Uyy and Uxy can be calculated by using the same matnx, where 

4 _ _ 2 

ai= 2.(^1 ”“^o) 

1=1 

4 

^2= 

1=1 

4 

= -yo) 

1=1 

4 __ _ 

Cl = 2(UXi -UxqXx, -Xq) 

1=1 

C2== Z(U3(i“Uyo)()^-yo) 

1=1 


Following these expressions any dependent vanable at the east, west, north and south 
faces of any cell for second order upwinding method can be wntten as 
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tej = <t)u + Uxij (xe- X .j)+ Uyy (ye- y u) 

<t>NlJ ~ ^t’lj ■** U xij (XN" X Ij)+ U yy (YN" Y y) 

(t>wy = <l)u + u xy (xw- X y)+ U jy (yw Y u) 
(J)Sy = (t)y + U xy(XS- X y)+ U yy(ys- Y Ij) 

where, 

Xe (X]H-i 

Xn =(x,+ 1j+1+X,j+i)/ 2 
Xw "(Xy ■^X,j+l)/2 

Xs ~(Xy ■^Xi+lj)/2 


The next step is to calculate the inviscid fluxes at the cell faces throughout the control 
volume The inviscid flux is approximated by the AUSM scheme as presented by Liou 
and Steffen (1993) For this scheme the normal Mach number at each face is calculated 
from the normal velocity vectors In each face of the cell this normal Mach number is 
used to calculate the normal flux For this the inviscid normal flux is divided into the 
convective and the pressure components given as 

F =[pu, puu , pvu, (E+p)u ]^+[0, p, 0, 0]^=F^‘^^+[0, p, 0, 0]^ 

G=[pu, pvu , pvv, (E+p)v ]’’+[0, 0, p, 0]’^=G^‘^+[0, 0, p, 0]’^ 

The effective way of writing the convective flux at any face of the cell is given as 
F ^'\/2 =M"i/ 2 [pa, pua, pva, (E-+p) a f l/r 
G^‘'\/ 2 =M"i/ 2 [pa, pua, pva, (E+p) a f ur 

where, 

(•)L/R=(*)LlfM",/2>0, 

= (•) R, Otherwise 

The advective Mach number M"i /2 is computed from the split Mach numbers M"\ and 
M"'r These split Mach numbers are denved according to the Van-Leer splitting as 
M^^ = ±1/4(M" ±l)^ if 1 M" |< 1 
= 1/2(M"±|M"|), otherwise 
The first order pressure splittmg can be expressed as, 
p^ =p(l±M")/2, if|M"i<l, 
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=p(M"±|M"|)/2M", otherwise 

Following these the flux formula at any cell face can be expressed as. 


M/2= 


'^pu '' 
puu + p 
pvu 

^(E + p)Uy 




1/2 



rpa 1 


fpa ^ 



'^pa ^ 


pua 


pua 


I 1 X ^ n » A 

pua 


pva 

+ 

pva 


--iMj/2 |Ai/2 

pva 


^(E + p)aj 

L 

JE + p)a; 

R. 


v(E + P)ay 


+ 


^0 ^ 

(pl +pr) 

0 

vO 


n/2“ 


'^pv ^ 



fpa ] 


'^pa ^ 



^pa '* 

pvu 

^ % jf XI 


pua 


pua 


^ B jr XT. 1 1 

pua 

pw + p 

= 2^1/2 


pva 

+ 

pva 


"2 ^^1/2 1^1/2 

pva 

^(E + p)vj 

1/2 


^(E + P)a> 

L 

^(E + p)a^ 

R- 


^(E+p)a^ 


+ 


ro 

0 

(pl +Pr) 
0 


(3 14) 


where, Aj/a { } = { }r- { }l Here { }r denotes the nght side cell of the face considered 
and { }l denotes the left side cell of that face 



Fig 3 1 Cell configuration in Finite Volume Method 
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To reduce the time of the computation, on the west and south faces of the cell (ij) the 
respective fluxes are copied from the fluxes of east and north face of the (i-lj) th and (ij- 
1) th cells respectively 

To limit the values of the gradients of the dependent vanables at the boundary cells of the 
flow domain the effect of the fourth cells ((i=4,in Fig 3 2) for the lower boundary and 
(i=2, m Fig 3 2) for the upper boundary) have not been taken in account For the inviscid 
flow the boundary condition is imposed by making the normal component of velocity 
zero on the free surfaces There is only pressure term that remains on those boundanes 
For the viscous flow simulations making the velocity vector zero at the solid walls 
imposes the boundary condition At the inlet all the dependent vanables are imposed for 
the supersonic flow and only the pressure term is extrapolated from the first inlet cell for 
the subsonic flow keeping all the other flow vanables imposed at the inlet Similarly at 
the outlet of the flow domain all the vanables are extrapolated for the supersonic flow and 
only one of the dependent vanables is explicitly imposed keepmg all other vanables 
extrapolated from the external boundary cell for the subsonic flow 

3 3 2 Calculation of Viscous Denvatives 

The viscous denvatives for the non-orthogonal mesh is evaluated following the Deiwart 
(1975), by transforming the vanables on local co-ordinate system concurrent to the mesh 
onentation. Fig 3 2 


Y 



‘ ► X 

Fig 3 2 Cell configuration for calculation of stresses 
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The transformation is given by, 

d\ dx 5ri ^ 

_^<yi 

dy d^dy dr] dy 


Here, ^ us a dummy-dependent variable and are the local coordinates of the non- 
orthogonal mesh The above transformation in matnx form becomes. 


4 >x 


"Hx 




i 

>* 

j 



(3 16) 


where, 


dx 


and^j 


dx 


respectively 1 he other symbols have their usual meanings 


Now, 


’^xTlx" 


'x^y^' 

-1 

_ 1 

■yti-y^' 



.XTiyr,: 

J 



where, J = 


Ixtiyp 


= (x^yT,-x^y^) 


Hence we have. 


<t>x 

_ 1 

■yti-y^' 

<t>^ 


J 

_x^-x^_ 



Therefore, 

d^ A<|)^Ay^-A(|)T,Ay^ 
dx Ax^Ay^ -Ax^Ay^ 

d^ (A({)^AXt,-A({)t,Ax^) 
dy (Ax^Ay^-Ax^^Ay^) 


For east face 




(3 17) 


(3 18) 
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A0^ =em,+j^j -0m, 

A0r| = +0m,j^.j ~0m,j_j — 0ni,+j j_|)/4 

For north face, 

A0m^ =(0m,+i^_,+0m,+j^j+, -0m,_i^j -0m,_j j+j)/4 
A0T1 =6mi,j+l 

In the preceding expressions, (|» denotes the velocity components u and v respectively, 
which are cell-centered vanables 0m denotes the center co-ordmates of the finite volume 
cell 

<t>=(u, v), 6m =(Xm, ym) 

If 0, j denotes the co-ordinate of one comer point of the hexahedral finite volume cell, m 
Fig 3 2, then 

©m,, =(9i,j +Qi+l,j +6i+l,j+l +9i,j+l)/4 

where, 0 = (x, y) 

This treatment of viscous denvative usually results in central differences and maintains 
second-order accuracy 

3 3.3 Time integration 

Using a four-stage Runge-Kutta scheme, the system of differential equations, Eqs 3 12b, 
are advanced in time It can be wntten for the nth time level as, 

U° =U" 

=U°-ai AtR(U®) 

=u‘^-a2AtR(U*) 

=U®-a3 AtR(U^) 

=U®-a4 AtR(U^) 
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Un+1 ^u4 


where, at the (m+l)st stage, 

R(U"''')-~Q(U"^) (3 19) 

with ai=I/4, a 2 = 1/3, a^=- 1/2, 04 - 1 0 

Based on the Founer stabiht> analysis of a one-dimensional model problem, the four- 
stage Runge-Kutta time stepping proposed by Jameson etal (1981) appears to be very 
attractive and is employed for the present work The CFL number restnction, based on the 
stability analysis of the model one-dimensional problem, is 2^2, which is relatively 
higher than most of the explicit schemes 


3 3 4 Acceleration Techniques 
Local Time Stepping 

I 

For faster steady state calculation, the locally varying maximum allowable time step is 
used The actual time step limit m the present work is expressed as 

At = CFL( - ) , where Ak is the limit due to the convective terms, Atj is the limit 

At(, + Atjj 

due to diffusive terms, and CFL is the Courant-Fnednehs-Lewy number In particular, 

. n 

Ate = 


Atn =• 


a' 


Kt(Yp/pPr) 


2-2 

+ s„ 


(3 20) 


where, and are as follows, 
X^ = q + c 
=qS^ + cS 


(3 21) 


where, and are the cell normal vector along each gnd line direction, q = (u, v) is 

the velocity vector and c is the speed of the sound and Kt is a constant that has been set to 
2 5 based on numencal experiments 
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Residual Smoothing 

By using implicit smoothing of the residuals, the stability range and robustness, along 
with the convergence to steady state of the basic time-steppmg scheme can be increased 
As suggested by Jameson and Baker (1983), the residual smoothing for three-dimensional 
flow, can be applied m the form, 

=R.j (3 22) 


where, R,,j = ■^Q(U) is the onginal residuals for the explicit time-stepping scheme and 


Rj j represents the residuals after the sequence of smoothing in ^ and ri directions with 

the coefficients and The residual smoothing is applied in the second and fourth 

time step of the four-stage Runge-Kutta integration In Eq 3 22, it is necessary to solve a 
sequence of tndiagonal equations for separate scalar vanables The use of constant 
coefficients m the implicit treatment has been proven to be beneficial in extending the 
CFL number to a considerable higher level than that in general used m the basic time 
stepping scheme without residual smoothing, which enhances the convergence rate quite 
significantly This has been quite saftsfactoiy even for highly stretched meshes, provided 
additional dissipative support such as enthalpy damping is mtroduced However, the use 
of enthalpy damping, which assumes constant total enthalpy throughout the field, 
precludes the solution of problems with heat transfer effects Using the above-mentioned 
can eliminate the need for enthalpy damping vanable coefficients P^ and P^^ that account 


for the variation in cell aspect ratio For the factorization of Eq 3 22, effective 
expressions for the coefficients p^ and p^ can be denved following the procedure of 


Martinelh (1987), as 

P4 


=:max{y[(- 

4 

= max{-^[( 
4 


CFL h 
CFL* 

CFL H 
CFL* ^ + 


s) 


(3 23) 


where, CFL is the Local Courant Number used m the computational scheme and CFL is 
the maximum allowable Local Courant Number based on stability analysis of the explicit 
Runge-Kutta scheme and <!>,, are defined as 


€> 4 = 1 +(“^) 


0 


CD ^I+( Jlf 

}ji 

where and are the quantities* given by the equation 3 21 and a has taken to be2/3 for 
the two-dimensional simulations The quantity e used in the above Eqs 3 22, is used as a 
limiter Dimensional analysis earned out found the value of s lying between 0 10 and 
0 25 For the present work, a value of 0 12 is used based on numencal expenmentation 
that assures stability along with enhancing the convergence For the present four-stage 
scheme. CFL = 1 0 and CFL* = 2 5 are used 


3.4 Closure 

The two-dimensional filtered compressible Navier-Stokes equations and the numencal 
algonthm in discretization of governing equations m the finite volume formulation are 
presented in this chapter The numencal procedure, along with the different acceleration 
techniques such as local time stepping and implicit residual smoothing has been discussed 
elaborately 
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CHAPTER -4 

RESULTS AND DISCUSSION 


4.1 Introduction 

A flov. solve! has been developed based on the finite volume formulation for 
compicssiblc flows The algonthm includes piecewise linear least square reconstruction, 
the Advectne I’pstrcam Splitting Method (AUSM) and a second-order time marching 
technique Ihe first section of this chapter deals with the computation of flows for 
different test cases both for supersonic and subsonic flow regimes, which illustrate the 
accuracy of the solver The flows over a 4% thick circular arc bump for different 
conditions are computed Another test case is the flow field developed due to shock- 
boundaT>’ layer interaction over a flat plate As already mentioned, the pnncipal objective 
of the work is to test the applicability of the flow solver to capture the unsteady flows by 
Large Eddy Simulation Thus the second part deals with the computation of the flow over 
a backward facing step by LES. Considering the time frame, the study is restricted to the 
two-dimensional analysis 

4.2 The Euler and Navier-Stokes analyses of flow over a 4% thick circular arc 
bump: 

As an illustrative problem, the supersonic and subsonic flows past a 4% thick circular arc 
bump in a channel are computed Both the Euler and Navier-Stokes analyses are 
presented, for the free-stream Mach numbers of M«:=l 4 and Moc^O 5 Predicted surface 
Mach numbers are compared with the previously computed results of Ni (1982) 

4 2 1 Gnd generation and boundary conditions 

The computational mesh for the mviscid calculation consists of 70x50 gnd points as 
shown in the Fig 4.1 The gnd is algebraically generated H gnd The corresponding mesh 
for the viscous calculation is shown in Fig 4 2 A highly refined stretched grid is used 
near the upper and lower boundanes for resolving the viscous layers A gnd sensitivity 
test IS conducted for choosing the number of gnd points need for the Navier-Stokes 
solution The surface Mach numbers for three different gnds of sizes 60x40, 80x60 and 
100x80 are plotted respectively (Fig 4 3) The surface Mach number obtamed from the 
refined gnd 100x80 is shown to be almost mdistmguishable from that with 80x60 gnds 
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Hence a gnd of 100x80 is chosen inside the flow domain for the Navier-Stokes 
calculations 

Now, the boundary conditions are descnbed At the inlet all the vanables are imposed for 
the supersonic flow For the subsonic inlet, the total pressure and the total temperature are 
specified and the axial component of the inlet velocity is computed as a part of the 
solution At the exit, all the vanables are extrapolated for the supersonic flow and the 
pressure is imposed in case of the subsonic flow For the Euler solution the normal 
component of velocity is made to zero at the wall On the wall, no-slip conditions are 
imposed for viscous calculation 

4 2 2 Results from the Euler and Navier-Stokes calculations 

It has been decided to examine the effect of the higher order discretization of the 
convective terms Fig 4 4, 4 5 and 4 6 illustrate the computed flows over a bump m a 
channel for supersonic inlet with different orders of spatial discretization The 
discretization process has been descnbed in chapter 3 It can be seen that two oblique 
shocks have formed at both comers of the bump The leading edge shock is reflected back 
by the upper wall into the region of expanding flow field The interaction between the 
reflected and trailing edge shocks near the trailing edge is captured by the solver, 
although the accuracy depends on the order of spatial discretization It is found that the 
second-order accurate scheme is less dissipative than the first-order The first-order 
AUSM is highly dissipative. It has been also seen that there is no significant difference 
between the solutions obtained from the third- and second-order of the spatial accuracy 
(Fig 4 5 and 4 6). 

Figure 4 7 demonstrate the companson of the surface Mach numbers obtained from the 
Euler solutions where convective terms have been discretised by the first-, second- and 
third-order upwind schemes The trend is similar as obtained earlier There is a significant 
improvement of the solution when second-order of spatial accuracy for the convective 
terms has been used However, the increase in accuracy is very small between the results 
obtamed from the second- and third-order of accurate solutions Thus Ihe second-order 
spatial discretization of the convective terms following AUSM will be adopted for further 
calculations 
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Figure 4 8 represents the surface Mach number distribution over the bump placed m a 
channel with the inlet Mach number of 1 4 The results obtained from the Euler solutions 
of second-order AUSM is compared with the results presented by Ni (1982) It is evident 
from the figure that the present calculation is in good agreement with that of Ni 
Moreover, the interaction of the reflected and trailing edge shock near the trailing edge of 
the bump is better resolved here. 

So far the results presented v, ere the solution of Euler equations Now, the Navier-Stokes 
computation i e the viscous calculations over the same bump arc are presented for the 
free stream Mach number of I 4 

Figure 4 9 represents the Mach contours for the viscous calculation over the bump for 
supersonic flow There is no appreciable difference between the Navier-Stokes and the 
corresponding Euler solutions (Fig 4 5) except near the trailing edge region of the bump 
This is due to the formation of the boundary layer and its interaction with the shock Here 
the shock has reflected from the top of the boundary layer. This phenomenon is absent m 
the Euler solution, where the shock will reflect from the wall Moreover, the flow may 
also separates due to the shock-boundary layer interaction The surface Mach number 
distnbutions are plotted in Fig 4 3 for the second-order Navier-Stokes solution The 
nature of the curve is same as that of the Euler solution Thus the present computation 
indicates that there exists very little difference between the surface pressure predicted by 
the Euler and the Navier-Stokes equations 

The skin-frnction coefficient over the bump is presented in Fig 4 10 It is found that a 
sudden increase m the skm-fhction coefficient occurs at the trailing edge The increase of 
the skin-fiiction occurs due to the acceleration of the flow over the bump, which causes 
the increase in velocity gradient at the wall At the trailing edge of the bump the flow 
meets with the reflected shock which results m sudden decrease of velocity close to the 
wall This is the result of the sudden drop of skm-fiiction at the trailing edge The flow 
has a tendency to separate there The similar trend has been reported by Parthasarathi and 
Kallmders (1995) 
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Figure 4 11 indicates the Mach contours obtained from the steady state Euler solution of 
flow over the bump placed inside a channel for subsonic flow (M«=0 5) It has been seen 
that the Mach contours are quite symmetric over the bump, which indicates the accuracy 
of the piesented computation The surface Mach number distribution has also become 
symmetric over the bump (Fig 4 12) Both of these results are similar to that of Ni 
(1982) I'he corresponding Mach contours over the bump obtained from the viscous 
calculation at M,-=0 5 is illustrated in Fig 4 13 Here the symmetnc nature of the Mach 
contours has lost due to the boundary layer growth as the flow past over the bump 
Howes er the nature of the surface Mach number distnbution remains almost similar to 
that of the Fuler solution as shown in Fig 4 14 

4 2 3 Conclusions 

'fhe computed results of the Mach contours and the surface pressure distnbutions both 
from the Euler and Navier-Stokes equations are presented for the flow over a circular arc 
bump placed inside a channel The calculations have been made both for subsonic and 
supersonic flow conditions Predicted results indicate the accuracy of the solver to resolve 
the flow pattern It should be noted that the second-order accurate AUSM scheme gives 
accuracy acceptable for engineenng applications 

4.3 The shock-boundary layer interaction over a flat plate 

ITiis section deals with the analysis of the boundary layer separation due to the 
impingement of a shock wave over a flat plate The results are compared with the 
available expenmental results 

4 3 1 Grid generation and the boundary conditions 

The grid used for the present study is shown in Fig 4 15 Here also an algebraically 
generated H gnd is chosen with high exponental gnd stretching near the wall A flow 
tangency condition is applied on the upper surface No-slip boundary conditions are 
imposed at the solid wall The total pressure and the total temperature are specified at the 
inlet section and the axial component of the inlet velocity is computed as a part of the 
solution At the outlet all the dependent variables are extrapolated. A weak shock is 
directed towards the boundary layer at an angle 32 58° with the direction parallel to the 
wall 
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4 3 2 Anaiy ses of th e shock-boundary layer interaction 

The schematic view of the shock-boundary layer interaction is shown in Fig 4 16 It can 
he seen that a shock is formed at the leading edge of the plate The shock creates an 
adverse pressure gradient causing the growth of the boundary layer there The incident 
shock IS directed at a certain angle parallel to the wall This shock reflects from the top of 
the boundar>’ layer and an adverse pressure gradient is created It also results in further 
giowth of boundary layer 'I'he adverse pressure gradient may be strong enough to cause 
the boundary layer separation This flow again reattaches the wall downstream of the 
shock as shown in the figure 

The numencal simulation of the shock-boundary layer interaction problem is performed 
with 70x56 gnd points The Mach contours are presented in Fig 4 17 for free-stream 
Mach number of 2 0 The incident shock and its reflection from the top of the boundary 
layer are captured properly by the AUSM algorithm The shock impingement creates an 
adverse pressure gradient, which increases the boundary layer thickness The increase in 
boundary layer thickness is clearly visible from flie Mach contours The boundary layer 
growth IS also same as that of diagram A leading edge shock similar to the schematic 
diagram is also visible from the figure 

The velocity vectors close to the wall for die flow are depicted m Fig 4 17 It is found 
that the adverse pressure gradient due to the incident shock not only thickens the 
boundary layer but also results m the separation of flow The separation, recirculation and 
the reattachment of flow are clearly illustrated by flie computed velocity field 

The skm-fiiction coefficient is plotted in Fig 418 The negative value of the skin-fhction 
coefficient further justifies that the separation of flow occurs where the shock meets with 
the boundaiy layer The result is venfied with the expenmental data by Hakkinen et al 
(1959) From this comparison, it is mferred that the AUSM under predicts the skin- 
fnction coefficient in the separation region Outside this region the computed skin-friction 
matches reasonably well with the expenmental data 
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4 3 3 C'onclusion'' 

I he anah*'!*' jiniK**’-' that the solver developed is accurate enough to resolve the 
complex How fu’ld ■ *1' Hig shock-boundary layer interactions and flow separation 

4.4 Two-dinu’nsi»i»<*l anahses of flovv over the backward facing step 
I he upsv ind s>. iieiiif *. m lesoKc flow structure well as it is based on the physics of flow 
\anuus upuiiKi ha\c been developed for LES of compressible flows An 

attempt is ni.u'u’ f<' '!5oulate the unsteady turbulent flow by using such an algonthm 
namcil as M SM li Ji.ts iscen decided to explore the possibility to apply the scheme for 
1 1*S I he tliiw o\i‘! a hovkuard facing step with 8 9 expansion ratio is chosen as a test 
case toi 1 i'S voinputatom However, the free-stream Mach number was very low 
M„-0 128 I his was chosen because of the fact that the expenmental data for the 
un.stead) Hows with uubuicnee quantities were available to us for this case 


4 4 1 ( ind gencialion and the boundary conditions 

The length along the stieamwise direction of the computational domain that is used for 
LE.S calculation is liioscn fiom a tnal steady state Navier-Stokes solution The separation 
region of flow has Ibund to be confined within an approximate length of 7H downstream 
of the step (where, H is the step height) A length of 16H along the streamwise direction 
after the step is considered to be sufficient for this LES calculation The height of the 
computational domain is taken as 9H that equal to the height of the test section A length 
of 811 i.s taken ahead of the step to match the computed boundary layer with that of 
experimental I’alue before the step The region is divided into algebraically generated 
non-uniform computational mesh of size 100x100 m the x and y directions The gnd is 
clustered near all the solid boundaries to capture steep gradients of the flow variables near 
the wall The computation domain and the grid are shown in Fig 4 19 The gnd is so 
chosen that the non-dimensional normal distance (y ) maintained of the order of unity 
near the wall A relatively high stretching of meshes near the step along the axial 
direction is chosen to predict the vortex structure generated inside the separation region 
A relatively coarse non-uniform mesh spacmg is used outside these regions 


33 



A ver> lovi. inlet Mach number turbulent flow (M«=0 128 and ReH=3 1x10^) is simulated 
m this problem At inlet as the boundary layer thickness is known the boundary layer 
profile IS specified considering 1/7 power law The static temperature is specified for free 
stream c ondition I he density is extrapolated at the inlet plane after each time step The 
static piessure is calculated from this density At the exit plane the static pressure is 
imposed and all the othei flou vanables are extrapolated At the upper and lower walls 
no-shp boundary conditions are imposed along with a vanishingly normal temperature 
giadient 1 he pressure at the solid surface is extrapolated from the intenor boundary cell 

4 4 2 Analy ses of flow over the backward facing step 

The soh er is first allowed to run for laminar flow with the same number of gnd points 
and above-mentioned boundary conditions The steady flow field generated is assumed as 
initial flow field for LES calculation Now, the SGS model is turned on for the rest of the 
simulation, which are continued upto 10 5^10^ iterations A CFL number of 1 0 is chosen 
for file entire simulation The lowest possible time step within the computational domain 
IS used for this simulation The time step is found to be very small, which needs about 10® 
Iterations to complete one eddy turn over time (Table 1) The very low free stream Mach 
number forces to fix time step to a small limit from stability entenon Thus, only 10 
samples have been collected for processing the results Thus the time-averaged results 
may not properly match with the available experimental data due to the small number of 
collected dataset The simulation is performed in a Pentium III machine with a clock 
speed of 800 MHz The computational time needed to complete 10® iterations is about 
1 30 hour> 


Table 1 


List of flow 

Vanables 

Values 

Free-stream 

0 128 

Mach number (M«) 


Average stream- 

44 2 m/sec 

wise velocity (Vavg) 


Step height (H) 

0 0127 m 
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List of flow 
\'anablcs 

Values 

Fddy turn-over 

7 9 X 10’^ sec 

time iH iiti 


I line step 

7 4 X 10'**sec 

si/eJ^Afl 



instantancitus sticanihiiCN and iickis are presented at different instants of time as 

shown in I ip 4 2(i (a-i) and 4 21 (a-d) ITic complex flow structure, formation and 
shedding of \’oitices and its dynamics are well resolved As the expansion ratio is very 
small there is no visible vortex formation on the upper wall The instantaneous flow 
pattern indicates that the flow field is unsteady 

The accurac) of results can be venfied by companng the time-averaged quantities with 
the experimental icsults Dnver and Seegmiller (1985). Time averaged streamwise and 
cioss-streamwise components of velocities are compared at 5 different locations as shown 
in Fig 4 22. The measurements at station 'a' (x/H=-2 0) give the initial value of the 
velocity fields. 'ITie station 'b' (x/H=I 5), station 'c' (x/H=4 0) provide the velocity fields 
at the beginning and middle of the recirculation zone The station 'd' indicates the section 
where the timc-avcraged recirculation zone reattaches In station 'e' (x/H=12 0) the 
velocity profiles are found to be fully developed It should be noticed that there exist 
discrepancies between the computed and the experimental results particularly for cross- 
stream components inside the recirculation zone (at x/H=l 5 and x/H=4 0) In the other 
sections the velocity profiles match more or less correctly with the expenmental results 
Similar trend was reported by P Sagaut et al (1994) However they concluded that these 
discrepancies could be minimized by using a djmamic subgnd scale model where the 
Smagonnsky constant automatically adjusts itself within the flow field 

Profiles of stream-wise (Urm/U) ^nd cross-streamwise (Vnns/U) turbulence intensities and 

the mean Reynolds stress (uV/U^) are plotted at different locafions (x/H=l 5, x/H=4 0, 
x/H=5 5 and x/H=12 0) as shown in Fig 4 23 It should be noted that the location 'a' and 
'b' are inside the separation region The section 'c' is just outside the recirculation zone 
and section 'd' is the region where the flow becomes folly developed All these turbulent 
quantities as shown in Fig 4 23 are averaged in both space (filtered) and time It can be 
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seen that the prc.scnt computation is able to capture the existence of a peak of turbulence 
intensities coirespondmg to the existence of a spatally evolving mixing layer issued from 
the separation legion although there are some discrepancies The anomalous peaks near 
the solid walls have been resulted due to the two-dimensional nature of the flow The 
two-dimensioiud simulation tan not predict vortex structure properly Similar trend was 
reporteil h\ F S.maut ct al (1994) Moreover, the small number of data sets used for time 
averaging of the How \ anabies enhances the error m the results 

An important parameter ot the backward facing step flow is the time averaged 
recnculation /one length ITiis length can be calculated from the skm-fhction coefficient 
The reattachment length is equal to the distance from the edge of the step to the pomt 
where the value of the coefficient changes from negative to positive The calculated value 
of this length is about 5.1 times the step height whereas the expenmental value is 6 26 
times the step height Thus, the simulation with the SGS model leads to the severe 
underpiediction of this quantity 

4 4 3 Con clusio ns 

From the above analyses it can be concluded that the Large-eddy simulation by using 
AUSM algontbm can be able to predict die turbulent flow over the backward facing step 
It has been seen that the Smagonnsky model is moderately accurate to calculate the 
turbulent stresses inside the flow domain However, for firm conclusion, it should be 
extended to thiee-dimensional simulation 
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Fig 4 1 C'ompulational mesh for the channel flow with a 4% thick circular arc bump for the 
liuler calculation 
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x/c 


Fig 4 2 Computational mesh for the channel flow with a 4% thick circular-arc bump for the 
Navier-Stokcs solution 





k/c 


Fig. 4.4 Mach contours for first order Euler solution of supersonic flow in a channel with a 
4% thick circular arc bump, M.=1.4. 



x/c 


Fig 4 5 Mach conu^urs for second order Euler solution of supersonic flow in a channel 
with a 4% thick ciicular arc bump, M«==l 4 
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1‘ig. 4.7 C’omparii»en uf surface Mach number distribution between the first, second 
and third order accurate results over the circular arc bump. 





Fig. 4.9 Mach contours for Navier-Stokes solution over the 4% thick circular arc bump 







Fig 4. 1 3 Mach contours of subsonic flow (Mco = 0 5) over a 4% thick circular arc bump 
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Fig. 4.14 Surface Mach number distnbution over a circular arc bump from Navier-Stokes 
solution with M« = 05 
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'ig 4 15 Computational mesh over the flat plate for shock-boundary layer interaction 
problem 




Fig 4 17 Mach contours for shock wave and boundary layer interaction 
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Fig 4 18 Companson of skm fhction coefficient between the experimental and the 
predicted value for the shock boundary layer interaction problem 
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Fig 4 i 9 Brtect of shock wave on the boundar 


5 






56 







(a) 



59 



(d) 


Fig 4 22 Vorticity plots at different instant oftrnie downstream of the 
backwaid facing step 
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Fig 4 23 Average streamwise and cross-streamwise components of velocity 
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CHAPTER -5 


CONCLUSIONS 


5.1 (“ondu.sions 

A two-dinicnsional flow solver is developed in physical plane using AUSM algonthm to split 
the com cclive fluxes A least square minimization procedure is adopted to define the fluxes 
at the cel! faces Higher order spatial discretization is made considermg the neighbourmg 
cells The time integration is performed by Runga-Kutta algonthm. The local tune steppmg 
and implicit residual smoothing are used to accelerate the steady state solution The scheme 
IS second order accurate m time and space For unsteady calculation a global minimum time 
step is used At present the calculation is limited to the two-dimensional simulation 

A number of test cases performed illustrate that the solver is reasonably accurate to simulate 
the complex flow pattern for a widely varying flow conditions from the subsomc to the 
supersonic flow regime Finally an attempt is made to simulate the turbulent flow structure 
over a backward facing step by LES using the flow solver A simplified SGS model as the 
Smagonnsky SGS model is used to model the subgnd scale turbulence. It has been found that 
the solver is able to resolve the unsteady vortex structures generated down-stream of the step 
There are some discrepancies between the time-averaged computed and the experimental 
results The reasons are for two-dimensional simulation and also for the low number of 
computed data sets used to generate time-averaged results Time limitation of an M Tech 
thesis would not allow the investigator to proceed more However, with computed results, it 
can be concluded that the present algonthm may be used to compute practical flows by LES 

5.2 Suggestions for Future Work 

The present Large Eddy Simulation has been performed over a backward facmg step with 
two-dimensional formulation The two-dimensional simulation can not predict the vortex 
structure properly The analyses can be extended to three-dimension so as to capture the 
proper vortex stretchmg generated downstream of the step As it has been found that the 
Smagonnsky model leads to high dissipation a dynamic model may be adopted to model the 



subgnd scale turbulence Finally the solver may be used to simulate very high Mach number 
flow s where giadient limiters should be adopted for stability restnctions 
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